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The Fermi gas at unitarity is a particularly interesting system of cold atoms, being dilute and 
strongly interacting at the same time. It can be studied non-perturbatively with Monte Carlo 
methods, like the recently developed worm algorithm. We discuss our implementation and tests 
I of this algorithm and suggest a modification that increases its efficiency by reducing autocorrela- 

^ tions. We then show how the worm algorithm can be applied to calculate the critical temperature 

Q of an imbalanced Fermi gas (unequal number of fermions in the two spin components). We fi- 

nally present some results obtained with the modified algorithm, in the balanced as well as in the 
imbalanced case. 
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1. Introduction 

Lattice methods are useful for studying strongly interacting theories in particle as well as in 
condensed matter physics. When strong interactions render a perturbative study impossible, lattice 
field theory can provide a useful tool for numerical calculations. The Fermi gas at unitarity is a 
prominent example for such a strongly interacting system [1]. 

Fermionic matter is ubiquitous in nature, from the electrons in metals and semiconductors or 
the neutrons in the inner crust of neutron stars, to gases of fermionic atoms, like ^^K or ^Li that 
can be created and studied under laboratory conditions. Due to Fermi-Dirac statistics, a dilute 
system of spin-polarised fermions exhibits no interactions and can be viewed as an ideal Fermi gas. 
However, interactions become crucial when we are dealing with fermions of several spin species. 
Low-energy interactions are characterised by the scattering length a. An especially intriguing case 
is the Fermi gas at divergent scattering length - the unitary regime, in which the gas is dilute (range 
of potential <sC inteiparticle distance) and strongly interacting (interparticle distance ^ scattering 
length) at the same time. One key feature of this regime is universality: since all information 
about interactions is contained in the scattering length, and this length scale is no longer present at 
unitarity, the gas exhibits universal behaviour that only depends on two dimensionful parameters, 
temperature and density. 

Due to a lack of an exact theoretical description several approximate and numerical approaches 
have been tried to study the Fermi gas in the unitarity limit. An accurate result for the critical tem- 
perature has been obtained with the recently developed Diagrammatic Determinant Monte Carlo 
(DDMC) algorithm [2]. In the following we will first introduce the model and the algorithm and 
then present our modifications and results. 

2. Fermi-Hubbard model at finite temperature 

The Fermi-Hubbard model is the simplest lattice model for two-particle scattering. Its Hamil- 
tonian in the grand canonical ensemble is given by 



where £k = ^ L]=i(1 ~ cosfcj) is the discrete dispersion relation, and (cka) the time-dependent 
fermionic creation (annihilation) operator. The model assumes non-relativistic fermions of two 
species labelled by a (which we will call "spin up" and "spin down") with equal particle mass m. 
For the present we also assume equal chemical potential pt^ = pt^ = pt for the two spin species. 
The attractive contact interaction is characterised by the coupling constant U <0. The limit of 
infinite scattering length corresponds toU = —7.914, in units where m = 1 /2. We work on a 3D 
simple cubic spatial lattice with sites, periodic boundary conditions and lattice spacing set to 
unity. The time direction remains continuous. The continuum limit of this model can be taken by 
extrapolation to vanishing filling factor. 

To study the finite temperature behaviour consider the grand canonical partition function in 



= //O + ^1 = £ (ek - iUa )4t,CkcT + [/ £ 4 
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Figure 1: Diagrammatic expansion of the partition function 



the imaginary time interaction picture, 



(2.2) 



where jS is the inverse temperature, T^- the imaginary time ordering operator and Hi{x) = e'^^°H\e 
Expanding Z in powers of H\ , 



oo . p 



p=Q 



e /^^o ncf(xy,T;-)cT(X;,T;)c|(xy,Ty)q(xy,T^-) 



, (2.3) 



generates a series of Feynman diagrams, as shown in figure 1. Since we are ukimately interested 
in thermal expectation values of operators and thermal averages are calculated using the expansion 
of the partition function, it would be convenient to use this expansion as a probability distribution 
to generate configurations for Monte Carlo sampling. However, each fermionic loop contributes a 
minus sign, with the consequence that the diagrams in the series have different signs. For our pur- 
pose we need to rewrite the series as a sum of positive terms only. This can be done by considering 
all diagrams of order p as one entity, which can be represented as the product of two matrix de- 
terminants (one p X p matrix for each spin component). The partition function can then be written 
as 

Z = det {Sp) det (Sp) (2.4) 

where Sp denotes a vertex configuration (the spacetime positions of all vertices) and the matrix 
entries are free finite-temperature single -particle propagators AJ^(5'p) = Gj^^(x/ — xy,T, — Ty). If 
the chemical potential is equal for spin up and spin down fermions (the balanced case) we have 
detA^ detA^^ = | det Ap, so that all terms in the series are positive [3]. 

The physical observable in the focus of our study is the order parameter for the phase transition 
to superfluidity. We introduce the pair creation and annihilation operators P(x,t) = Cx|(t)cx|(t) 
and P^(x',t') = c^/|(f')Cx'j.('^')- the critical point the correlation function 



G2(xt;x't') = (T,P(x,t)P^(x',t'^ 



1 



Tr[T,P(x,T)P^(x',TO^"'^''] 



(2.5) 



is proportional to |x — x'| as |x — x'| — > oo, where T] 0.038 is the anomalous dimension for 

the U(l) universality class. The integrated correlation function 



Tti Jo Jo 



(2.6) 
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will be independent of the lattice size L at jSc = 1/7"^. This property can be used to determine the 
critical temperature. 

3. Worm algorithm 

The configuration space is sampled via a Monte Carlo Markov chain process: in each step 
one of the possible updates to another vertex configuration is proposed, and accepted with some 
probability given by the detailed balance equations. The requirements of detailed balance and er- 
godicity ensure that the produced configurations are indeed distributed according to the correct 
thermal probability distribution p{Sp) = ^(— ?7)''| det A(5p) given by the expansion of the parti- 
tion function. 

The diagrammatic expansion of Tr[TTP(x, t)P^ (x', x')e^^^] is similar to that of Z, but contains 
an additional pair of 2-point vertices at (x, t) and (x', t'). It is thus of advantage to sample these 
two series in the same simulation. In addition to sampling the regular 4-point diagrams we allow 
updates that insert the pair of 2-point vertices ("worm vertices") into the configuration space. This 
setup has several advantages. Firstly, the Monte Carlo estimator for the order parameter becomes 
very simple: it is the ratio of configurations with and without worm vertices. Secondly, all updates 
can now be performed through the worm vertices P and , which simplifies the setup. 

A detailed description of the individual updates can be found in the appendix of [2]. Here we 
will only give a brief summary. 

• Updates only concerning the worm vertices: 

- Worm creation/annihilation: insert/remove the pair P(x, t), P^(x',t') into/from the 
configuration. 

- Worm shift: shift the P^(x', t') vertex to other coordinates. 

• Updates of the regular 4-point vertices: adding/removing a 4-point vertex (changes the dia- 
gram order). 

- Diagonal version: add or remove a random vertex. 

- Alternative using worm: move the P{x, x) vertex to another position and insert a 4-point 
vertex at its old position. The new coordinates are chosen in a way that tends to prolong 
existing vertex chains. In this case an update can only happen when the pair of 2-point 
vertices is present. 

The worm setup, as proposed in [2], leads to much higher acceptance rates than the regular 
"diagonal setup". The idea behind the worm setup is that at low densities the major contribution 
comes from multi-ladder diagrams, these are configurations where the vertices are arranged into 
several vertex chains. The typical size of a chain depends on the parameters of the system. To 
favour the creation of vertex chains the worm update uses the 2-point vertex P to add (or remove) 
4-point vertices in a small spacetime region, while P gets shifted to new coordinates each time. 
At low densities the new coordinates of P will be chosen according to a probability distribution 
that favours zero or small spatial shifts and small temporal shifts in the direction of positive T. 
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The removal update always attempts to remove the nearest neighbour of P, which means that due 
to detailed balance the addition update can only be accepted if the added vertex is the nearest 
neighbour of the shifted P-vertex. This nearest neighbour condition is crucial for achieving high 
acceptance rates. 

4. Modifications of the algorithm 

In our study we found that although the worm type addition and removal updates have high 
acceptance rates, it is at the cost of strong autocorrelations. It is most efficient for the algorithm 
to successively add and remove the same vertices, so that the configuration does not change sig- 
nificantly, even after many successful updates. To illustrate this compare the measurements of the 
interaction energy (which is proportional to the diagram order) in the diagonal and the worm setup 
(figures 2 and 3). Both simulations used the same parameters and a comparable number of MC 
steps. 

Because of the large measurement en^or due to autocorrelations the worm setup is effectively 
less efficient than the standard diagonal setup. These new insights make a modification of the 
worm algorithm necessary. The goal is hereby to combine the advantages of the diagonal setup 
(weak autocorrelations) with the ones of the worm setup (high acceptance rates). To achieve this 
we propose a second type of addition/removal updates: 

• Choose a random 4-point vertex from the configuration (which will act as a worm for this 
step). 

• Addition: add another 4-point vertex on the same lattice site and in some time interval around 
the worm. 

• Removal: remove the nearest neighbour of the worm vertex (implies that addition can only 
be accepted if the new vertex is the nearest neighbour of the worm). 

This setup still prolongs existing vertex chains, but autocon^elations are significantly reduced since 
the worm changes with every update. This new type of updates can of course only be employed 
in addition to the regular diagonal addition and removal updates. It works regardless if the pair 
of 2-point vertices is present or not in the configuration (the worm addition/removal updates can 
only take place when the 2-point vertices are present). The acceptance rates for this update are 
comparable with those for the regular worm updates. 

5. Imbalanced Fermi gas 

The DDMC algorithm presented in the previous sections relies strongly on the assumption of 
equal densities of the two fermion species. This assumption allows us to write the partition function 
(2.4) as a sum of positive terms only, and consequently to use it as a probability distribution for 
Monte Carlo sampling. We now present a generalisation of the algorithm to the imbalanced unitary 
Fermi gas (/i| ^ pt^). The imbalanced case is especially interesting for a variety of reasons. A 
much richer structure of the phase diagram can be observed in this case. The superfluid state 
has been found to be remarkably stable, however it is also known that at some critical imbalance 
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(a) Diagonal setup (b) Worm setup 

Figure 2: The first 100000 measurements of the interaction energy (a measurement takes place every 100 
MC steps). Strong autocorrelations are visible in the worm setup. 




(a) Diagonal setup (b) Worm setup 

Figure 3: The blocking error analysis of the interaction energy (lines to guide the eye). The blocked error is 
much higher in the worm setup and continues increasing even for large block sizes. 



superfluidity must break down completely [4]. Several experimental studies are already available 
[5], but to our knowledge numerical studies have so far been only performed at zero temperature 
[6]. 

Our goal is to study how an imbalance will affect the critical temperature of the system. In this 
case a sign problem is present: the function p{Sp) = |(— ?7)^detAl^(5'p) detA^(5'p) is no longer 
positive for all configurations Sp and can thus not be used as a probability distribution. Several 
methods of dealing with sign problems of this kind are known from lattice QCD, where the intro- 
duction of a chemical potential renders the fermionic determinant complex. The most straightfor- 
ward one is the "phase quenched method", which reduces to a "sign quenched method" in our case. 
We can write p{Sp) = |p(5p)|sign(5p) and use the positive function |p(5'p)| as the new probability 
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distribution. For a thermal average this means 

^ ZnSp)p{Sp) ^ ZnSp)\p{Sp)\^ign{Sp) ^ (Xsign)ipi 
^ LPiSp) L\p{Sp)\^ign{Sp) (sign)ipi 

This representation of a thermal average in terms of the new probability distribution |p(5p)| is 
mathematically equivalent to the usual thermal average. However, numerical errors can become 
very large if (sign)|p| « 0, as it happens for the expectation value of the phase in QCD. Our studies 
have shown that for the unitary Fermi gas the sign remains very close to unity for a wide range of 
imbalances, so that sign quenching is applicable and accurate values for the critical temperature 
can be obtained. Some preliminary data for the balanced and imbalanced case will be presented in 
the next section. 



6. Results 

With the modified worm algorithm we were able to reproduce several values of Tc at different 
filling factors v for the balanced case. Figure 4 shows a typical order parameter analysis. Since the 
integrated correlation function R{L,T) defined in equation (2.6) becomes lattice size independent 
at the critical point, the R{L, T) curves for different L are expected to cross at j8c = 1 /T^. Due to 
finite size effects the curves do not cross in exactly one point, but renormalisation group analysis 
can be applied to extrapolate to the infinite volume limit [2]. For the point in 4(a) we obtain 
vV3 = 0.542 ±0.003 and Tc = (0.087 ± 0.002)£'f, where Ep is the Fermi energy. This agrees 
well with the results presented in [2]. In figure 4(b) we present the crossing of the order parameter 
lines for the unitary Fermi gas with an imbalance of |A/x| = (0.0398 ±0.0007)£'/7. In this case the 
expectation value of the sign was found to be between 0.98 and approximately 1 for lattice sizes 




(a) balanced (b) imbalanced 

Figure 4: The integrated correlation function R{L, T) is plotted versus the inverse temperature (5 = 1/7 for 
different lattice sizes. The lines cross near the critical point, where the order parameter becomes lattice size 
independent. 
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L = 12,8,6. We can see that the overall errors are sufficiently small. For the point in 4(b) we obtain 
= 0.512 ±0.005 and = (0.084 ±0.005)£f. 

Our studies so far indicate that the sign quenched method can be successfully applied to imbal- 
ances of up to about 0.3Ef. For very large imbalances the sign average becomes close to zero, so 
that reliable results can no longer be obtained. Further measurements for a wider range of densities 
and imbalances are in progress. 

Acknowledgments 

This work has made use of the resources provided by the Cambridge High Performance Com- 
puting Facility. OG is supported by the German Academic Exchange Service (DAAD), the Engi- 
neering and Physical Sciences Research Council (EPSRC) and the Cambridge European Trust. 

References 

[1] for a review see e.g. S. Giorgini, L. P. Pitaevskii and S. Stringari (2007), [arXiv : 07 6 . 33 60] 
[2] E. Burovski, N. Prokof 'ev, B. Svistunov and M. Troyer, New J. Phys. 8, 153 (2006) 

[3] A. N. Rubtsov, V. V. Savkin and A. I. Lichtenstein, Phys. Rev. B 72, 035122 (2005), 

[cond-mat/041134 4] 

[4] D. T. Son and M. A. Stephanov, Phys. Rev. A 74, 013614 (2006) 

[5] M. W. Zwierlein, A. Schirotzek, C. H. Schunck, W. Ketterle, Science 311, 492 (2006); 
G. B. Partridge, W. Li, R. I. Kamar, Y. A. Liao, R. G. Hulet, Science 311, 503 (2006); 
Y. Shin, M. W. Zwierlein, C. H. Schunck, A. Schirotzek, W. Ketterle, Phys. Rev Lett. 97, 030401 
(2006); 

G. B. Partridge, W. Li, Y. A. Liao, R. G. Hulet, M. Haque and H. T. C. Stoof, Phys. Rev Lett. 97, 
190407 (2006) 

[6] J. Carlson and S. Reddy, Phys. Rev Lett. 95, 060401 (2005), [cond-mat/050325 6]; 

L. M. Jensen, J. Kinnunen and P Torma, Phys. Rev A 76, 033620 (2007), [cond-mat/0604424] 



8 



